0.1 Packages

library(tidyverse)
library(broom)

1.1 Life Distributions

All technologies, operations, etc. have a ‘lifetime distribution’. If you took a sample of, say, cars in New York, you could measure how long each car functioned properly (its life-span), and build a Lifetime Distribution from that vector.

# Let's imagine a normally distributed lifespan for these cars...
lifespan <- rnorm(100, mean = 5, sd = 1)

The lifetime distribution is the probability density function telling us how frequently each potential lifespan is expected to occur.

# We can build ourself the PDF of our lifetime distribution here
dlife <- lifespan %>% density() %>% tidy() %>% approxfun()

In contrast, the Cumulative Distribution Function (CDF) for a lifetime distribution tells us, for any time \(t\), the probability that a car will fail by time \(t\).

# And we can build the CDF here
plife <- lifespan %>% density() %>% tidy() %>% 
  mutate(y = cumsum(y) / sum(y)) %>% approxfun()

Having built these functions for our cars, we can generate the probability (PDF) and cumulative probability (CDF) of failure across our observed vector of car lifespans, from ~2.7 to ~7.44.

Reliability or Survival Analysis is concerned with the probability that a unit (our car) will still be operating by a specific time \(t\), representing the percentage of all cars that will survive to that point in time. So let’s also calculate 1 - cumulative probability of failure.

mycars <- tibble(
  time = seq(min(lifespan), max(lifespan), by = 0.1),
  # Get probability of failing at time time
  prob = time %>% dlife(),
  # Get probability of failing at or before time t
  prob_cumulative = time %>% plife(),
  # Get probability of surving past time t 
  # (NOT failing at or before time t) 
  prob_survival = 1 - time %>% plife())

Let’s plot our three curves!

ggplot() +
  # Make one area plot for Cumulative Probability (CDF)
  geom_area(data = mycars, 
            mapping = aes(x = time, y = prob_cumulative, 
                          fill = "Cumulative Probability"), alpha = 0.5) +
  # Make one area plot for Relibability
  geom_area(data = mycars, 
            mapping = aes(x = time, y = prob_survival, 
                          fill = "Reliability (Survival)"), alpha = 0.5) +
  # Make one area plot for Probability (PDF)
  geom_area(data = mycars,
            mapping = aes(x = time, y = prob, 
                          fill = "Probability"), alpha = 0.5) +
  theme_classic() +
  theme(legend.position = "bottom") +
  labs(x = "Lifespan of Car", y = "Probability",
       subtitle = "Example Life Distributions")

Figure 2. Life Distributions of a Fleet of Cars

This new reliability function allows us to calculate 2 quantities of interest:

• expected (average) number of cars that fail up to time \(t\).

• total cars expected to still operate after time \(t\).

1.2 Example: Airplane Propeller Failure!

Suppose Lockheed Martin purchases 800 new airplane propellers. When asked about the failure rate, the seller reports that every 1500 days, 2 of these propellers are expected to break. Using this, we can calculate \(m\), the mean time to fail!

\[ \lambda \ t = t_{days} \times \frac{2 \ units}{1500 \ days} \ \ \ and \ \ \ m = \frac{1500}{2} = 750 \ days \] This lets us generate the failure rate \(F(t)\), also known as the Cumulative Distribution Function, and we can write it up like this.

\[ CDF(t) = F(t) = 1 - e^{-(2t/1500)} = 1 - e^{-t/750} \] What’s pretty cool is, we can tell R to make a matching function fplane(), using the function() command.

# For any value `t` we supply, do the following to that `t` value.
fplane = function(t){ 1 - exp( -1*(t / 750)) }

Let’s use our function fplane() to answer our Lockheed engineers’ questions.

1. What’s the probability a propeller will fail by t = 600 days? By t = 5000 days?

fplane(t = c(600, 5000))

## [1] 0.5506710 0.9987274

Looks like 55% will fail by 600 days, and 99% fail by 5000 days.

2. What’s the probability a propeller will fail between 600 and 5000 days?

fplane(t = 5000) - fplane(t = 600)

## [1] 0.4480563

~45% more will fail between this period.

3. What percentage of new propellers will work more than 6000 days?

1 - fplane(t = 6000)

## [1] 0.0003354626

0.03% will survive past 6000 days.

4. If Lockheed uses 300 propellers, how many will fail in 1 year? In 3 years?

# Given a sample of 300 propellers,
n <- 300
# We project n * fplane(t = 362.25) will fail in 1 year (365.25 days)
# that's ~115 propellers.
n*fplane(t = 365.25)

## [1] 115.6599

# We also prject that n * fplane(t = 3 * 362.25) will fail in 3 years
# that's ~230 propellers!
n*fplane(t = 3*365.25)

## [1] 230.3988

Pretty powerful!

2.1 Multiplication Rule

• probability that \(n\) units with a reliability function \(R(t)\) survive past time \(t\) is multiplied, because of conditional probability, to equal \(R(t)^{n}\)..

Figure 4. Oops.

For example, there’s a 50% chance that 1 coffee cup breaks at local coffeeshop Coffee Please! every 6 months (180 days). Thus, the mean number of days to cup failure is \(m = \frac{180 \ days}{ 1 \ cup \times 0.50 \ chance} = 360 \ days\), while the relative frequency (probability) that a cup will break is $ = $.

fcup = function(days){ 1 - exp( -1*(days/360))}
# So the probability that 1 breaks within 100 days is XX percent
fcup(100)

## [1] 0.2425349

And let’s write out a reliability function too, based on our function for the failure function.

# Notice how we can reference an earlier function fcup in our later function? Always have to define functions in order.
rcup = function(days){ 1 - fcup(days) }
# So the probability that 1 *doesn't* break within 100 days is XX perecent
rcup(100)

## [1] 0.7574651

But the probability that two break within 100 days is…

fcup(100) * fcup(100)

## [1] 0.05882316

And the probability that 5 break within 100 days is… very small!

fcup(100)^5

## [1] 0.0008392106

2.2 Compliment Rule

• The probability that at least 1 of \(n\) units fails by time \(t\) is \(1 - R(t)^{n}\).

So, if Coffee Please! buys 2 new cups for their store, the probability that at least 1 unit breaks within a year is…

1 - rcup(days = 365.25)^2

## [1] 0.868555

While if they buy 5 new cups for their store, the chance at least 1 cup breaks within a year is…

1 - rcup(days = 365.25)^5

## [1] 0.9937359

3.1 Hazard Rate Function

But if a unit has survived up until now, shouldn’t its odds of failing change? We can express this as:

\[ P(Fail \ Tomorrow | Survive \ until \ Today) = \frac{ F(days + \Delta days) - F(days) }{ \Delta days \times R(days)} = \frac{ F(days + 1 ) - F(days) }{ 1 \times R(days)} \] Local coffeeshop Coffee Please! also has a lucky mug, which has stayed intact for 5 years, despite being dropped numerous times by patrons. Coffee Please!’s failure rate suggests they had a 99.3% chance of it breaking to date.

# we call this the Hazard Rate Z
zcup = function(days, plus = 1){ 
  ( fcup(days + plus) - fcup(days)  )  / (plus * rcup(days) )
}

3.2 Accumulative Hazard Function

• \(H(t)\): total accumulated risk of experiencing the event of interest that has been gained by progressing from 0 to time \(t\).

• the (instantaneous) hazard rate \(h(t)\) can grow or shrink over time, but the cumulative hazard rate only increases or stays constant.

hcup = function(days){ -1*log( rcup(days) ) }
# This captures the accumulative probability of a hazard (failure) occurring given the number of days past.
hcup(100)

## [1] 0.2777778

3.3 Average Failure Rate

The hazard rate \(z(t)\) varies over time, so let’s generate a single statistic to summarize the distribution of hazard rates that \(z(t)\) can provide us between times \(t_{a} \to t_{z}\). We’ll call this the Average Failure Rate \(AFR(T)\).

afrcup = function(t1,t2){ (hcup(t2) - hcup(t1) ) / (t2 - t1)}
afrcup(0, 5)

## [1] 0.002777778

# Or write it as....
afrcup = function(t1,t2){ (log(rcup(t1)) - log(rcup(t2)) ) / (t2 - t1)}
afrcup(0, 5)

## [1] 0.002777778

# And if we're going from 0 to time t,
# it simplifies to...
afrcup = function(days){ hcup(days) / days }
afrcup(5)

## [1] 0.002777778

# or to this
afrcup = function(days){ -1*log(rcup(days)) / days }
afrcup(5)

## [1] 0.002777778

When the probability for a time \(t\) is less than 0.10, \(AFR = F(t) / T\). This means that \(F(t) = 1 - e^{-T \times AFR(T)} \approx T \times AFR(T) \ \ when \ F(T) < 0.10\).

afrcup = function(days){ fcup(days) / days }
afrcup(5)

## [1] 0.002758577

# and this is approximately....

3.4 Table of Functions

4.1 Failure Functions

First, let’s write the failure function f(t).

# Write failure rate
f = function(t){ 1 - exp(-(t/2000)^0.5) }

Second, let’s write the hazard rate z(t), for a 1 unit change in t.

# Write hazard rate
z = function(t, change = 1){ 
  # Often I like to break up my functions into multiple lines;
  # it makes it much clearer to me.
  # To help, we can make 'temporary' objects; 
  # they only exist within the function.
  
  # Get change in failure function
  deltaf <- (f(t+change) - f(t) ) / change
  # Get reliability function
  r <- 1 - f(t)
  # Get hazard rate
  deltaf / r
}

Third, let’s write the average hazard rate afr(t1,t2).

afr = function(t1,t2){
  
  # Let's get the survival rate r(t)
  r1 = 1 - f(t1)
  r2 = 1 - f(t2)
  
  # Let's get the accumulative hazard rate h(t)
  h1 = -log(r1)
  h2 = -log(r2)
  
  # And let's calculate the averge failure rate!
  afr = (h2 - h1) / (t2 - t1)
  
  # And return it!
  afr
}

4.2 Conversion Functions

Fourth, let’s write some functions to convert our results into %/K and PPM/K, so we can be lazy! We’ll call our functions pk() and ppmk().

# % per 1000 hours
pk = function(rate){ rate * 100 * 10^3 }

# PPM/1000 hours
ppmk = function(rate){ rate * 10^9 }

4.3 Converting Estimates

Let’s compare our hazard rates when t = 10, per hour, in % per 1000 hours, and in PPM per 1000 hours.

# Per hour... Ew. Not very readable.
z(t = 10)

## [1] 0.003445358

# % per 1000 hours.... Wheee! Much more legible
z(t = 10) %>% pk()

## [1] 344.5358

# PPM per 1000 hours.... Who! Big numbers!
z(t = 10) %>% ppmk()

## [1] 3445358

Finally, let’s calculate the Average Failure Rate between 1000 and 10000 hours, in %/K.

# Tada! Average Failure Rate from 1000 to 10000 hours, in % of units per 1000 hours
afr(1000, 10000) %>% pk()

## [1] 16.98846

# And in ppmk!
afr(1000, 10000) %>% ppmk()

## [1] 169884.6